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Abstract 

We develop a method of stochastic differential equation to simulate electron acceleration at astrophysical shocks. Our 
method is based on Ito’s stochastic differential equations coupled with a particle splitting, employing a skew Brownian 
motion where an asymmetric shock crossing probability is considered. Using this code, we perform simulations of electron 
acceleration at stationary plane parallel shock with various parameter sets, and studied how the cutoff shape, which is 
characterized by cutoff shape parameter a, changes with the momentum dependence of the diffusion coefficient /3. In the 
age-limited cases, we reproduce previous results of other authors, a ~ 2/3. In the cooling-limited cases, the analytical 
expectation a « /3 + 1 is roughly reproduced although we recognize deviations to some extent. In the case of escape- 
limited acceleration, numerical result fits analytical stationary solution well, but deviates from the previous asymptotic 
analytical formula a ~ /3. 
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1. Introduction 

Mechanism of particle ac celeration is sti l l unknown . 
Diffusive shock acceleration ( Krvmski I Il977t iBell Il978l: 
Blandford fe Ostrikeil 1978 1 is the most plausible if strong 
shock waves exist as in young supernova remnants (SNRs). 
We have not yet well constrained model parameters, namely 
magnetic field strength and degree of magnetohydrody¬ 
namic turbulence, although there are obser vat ional claims 
of turbulent, amplif i ed field in you n g SNRs (Vink fe LaminsJ . 
200g|^ Bambajd ak^ 20031 l2005al lb: Yamazaki et al. . 2004 : 


field strength. 


IJchivarna et al. . 2007tl . These are important to estimate 
maximum attainable energy of both electrons and nuclei 


fe.g.. lYoshida fe Yanagital |1997|) . Yamazaki et al. (|2013lf 


proposed that cutoff shape of electron spectrum around 
the maximum energy E max may provide us important in¬ 
formation on the cosmic-ray acceleration at young SNRs. 
They related the cutoff shape parameter a, which is defined 
by N(E) oc exp[— (U/U max ) a ], to the energy dependence 
of the electron diffusion coefficient /3 (that is, K oc E@) 
in each case where the maximum electron energy is de¬ 
termined by SNR age, synchrotron cooling and escape 
from the shock. They found that if the power-law in¬ 
dex of the electron spectrum is independently determined 
by other observations, then the cutoff shape parameter 
can be constrained by near future hard X-ray observations 
such as Nuclear Spectroscopic Telescope A rray (NuSTAR) 


(Hai lev et all l201flt lHarrison et all 1201,311 and ASTRO-H 


( Takahashi et al. , 201f)h and/or CTA ( Actis et all 1201 lh . 


These X-ray and gamma-ray observations will be impor¬ 
tant for the estimate of /3 as well as E max and the magnetic 


In analysis of Yamazaki et al. ( 2013h 


they assumed 

relations between a and /3 as a = 2/3, /3 + 1 and /3 in 
the case of age-limited, cooling-limited and escape-limited 
acceleration, respectively. The formula a = 2/3 in the 
age-limited case has been based on numer ical simulation 
( Kato fe Takaharal 120031 : iKang et al~ . 20091) . while the oth¬ 


ers are obtained analytically on the assumption of station¬ 
ary state, and they are not yet confirmed numerically. In 
this paper, we study the cutoff shape of the electron spec¬ 
trum by numerically solving the transport equation de¬ 
scribing diffusive shock acceleration, and study whether 
the above relations are right or not. 

We use a numerical method for solving cosmic-ray trans¬ 
port equation (so-called, diffusion-convection equation), 
which was proposed by Achterberg fe Kriillsl ( 19921) . This 


method is based on the equivalence between the Fokker- 
Planck equation and the Ito stochastic differential equa- 


lowed for various situations (Kriills & Achterberg. 

1994; 

Yoshida & Yanagita. 1994: Marcowith & Kirkj. 19991: 

Marcowith & 

2010/ Schure et alj, 2010|). It should be noted that the 


SDE method has an advantage if the transport equation 
has to be solved in multi-dimensions. In practice, the im¬ 
portance of upstream inhomogeneity for understanding of 
cosmic-ray acceleration at supernova rem na nts h as been 
pointed out by various authors (e.g., Ilnoue et~all l2012h . 
In this case, it is clear to consider the particle acceleration 
in three dimensions. 

The simple-minded application of the SDE method has 
problems in actual numerical integration. First, (5-functions 
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appear in SDE if we apply it to the shock front, where the 
background fluid velocity as well as the diffusion coeffi¬ 
cient have a sudden jump. In ord er to avoid this, the shock 
structure is artificially smoothed ( Achterberg fe Kriills . 1992h . 

However, even in this case, the time step has to be small 
enough for the simulated particles not to miss the sharp 
gradient at the shock front, which significantly slows down 
the simulation. Furthermore, in actual simulation time, 
approximation of the smooth shock transition cause s incor¬ 
rect p article spectrum. This difficulty was s olved by Zhanel . ._. 

j200d) who used the skew Brownian motion ( Harrison fc Shepnl f c ^ ua '^ on ^ 
Il98lh which can be solved by a scaling method that elimi- gp 


and c are Thomson cross section, 
velocity of light, respectively. 
Introducing new quantities, 


P 


mass of electron and 


( 4 ) 


and 


F(x, u, t) = p 3 f{x,p,t) 


becomes the Fokker-Planck form, 


( 5 ) 


nated the (5-functions in the SDE. Other numerical schemes 


to resolve this problem have also b een proposed (IMarcowith fe Kirk , 
Il999t lAchterberg fe Schurel . 201 ill . Second problem is that 
a large dynamic range in particle momentum causes low 
statistical accuracy at large momenta. This difficulty was 
also res ol ved by e mpl oying a particle splitting technique 
( Yoshida fe Yanagital . 1 1 99411 . 

In this paper, we first attempt to perform simul ations 
of elec tron acceleration incorporating both methods of Zhanel 
( 200d l and particle splitting. Owing to newly developed 
code, simulated spectra have cutoff shape accurate enough 
to be compared with analytical formulae. As a first step, 
we focus on the cases of one-dinrensional plane shock. Ex¬ 
tended studies for more complicated cases such as time de¬ 
pendent free escape boundary, nonuniform magnetic fields, 
and/or multi-dimensional systems (including spherical shock 
geometry) are simple but remained as future works. 


dt 


d_ 

dx 


d_ 

du 


<)K 
dx 
1 dv 
3 dx 


dx 2 


(KF) 


+ Psynl F 


= 0 


( 6 ) 


This equation is equivalent to the following SDEs of the 
Ito form: 


dx = 


du = — 


dK 

dx 


-- 1 dt + V2KdW 


1 dv 
3 dx 


fisyn'y J dt 


( 7 ) 


( 8 ) 


where dW is a Wiener process given by the Gaussian dis¬ 
tribution: 


P(dW) = 


1 


V^ndt 


exp(— dW 2 /2dt) 


( 9 ) 


2. Basic Equations and Numerical Method 

2.1. Basic equations 

In this paper, we consider one-dimensional system, that 
is, all quantities depend on the spatial coordinate x. The 
diffusion-convection equation with energy-loss process is 
given by 


d£ 

dt 


JLf 

dx A dx 


p 


l_d_ 
2 dp 


pdv dp i 2 
~3 


= 0 


(1) 


where f(x,p,t) is the distribution function for electrons, 
and p is the electron momentum. Functions v(x) and 
K(x,p) are background velocity field and the spatial diffu¬ 
sion coefficient of the electrons, respectively. In this paper, 
we consider the synchrotron cooling. Then, the loss term 


becomes 



dp 

dt 

-psynlP , 

( 2 ) 

where 



fisyn — 

cr^B 2 

67 mi e c 

( 3 ) 


and 7 = ( p/m e c ) 2 + 1 is the electrons’ Lorentz factor, 
and B is the magnetic field. Physical constants, ctt, m e 


Numerical simulation by SDEs is much faster than that 
with the usual Monte Carlo method and is much easier 
than solving the original Fokker-Planck equation, because 
the SDEs are ordinary differential equations. 

2.2. Method of Zhang (2000) 

The application of the SDEs, equations 0 and (ED, 
for the study of electron acceleration at the shock is not 
simple, because the velocity field v(x) has a sudden jump 
at the shock front, so that dv/dx in equation ED contains 
5-function. Similarly, if the diffusion coefficient also be¬ 
haves discontinuously at the shock front, then dK/dx in 
equation 0 also contains the 5-function. We take the co¬ 
moving frame with the shock which is located at x = 0 
a nd w e define x < 0 as upstream region. Following I Zhang 
iteoooli . we decompose the velocity field v and the diffusion 
coefficient K into two parts: 


, , , , AH . 

v{x) = v c {x) + sign (a:) , 

( 10 ) 

K(x) = K c (x) + 4^-sign(:r) , 

( 11 ) 
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where A.V = u(0 + ) — v(0~) and A I\ = AT(0 + ) — K{0~), 
and sign(a;) is the sign of x. Functions v c (x) and K c (x) 
are continuous for arbitrary x (including x = 0). We scale 













































the x coordinate accordin g to its sign in the following way 
( Harrison fe Sheml . Il98lh : 


y = xs(x) = x x 


a 

1 

2 


(x < 0) 
(x = 0) 


x 


( 12 ) 


1 — a (x > 0) 


0 (xiX i+ 1 > 0) 

x i+ i/(a - 1) ( xi > 0, x i+ i < 0) 

Xi+i/a ( Xi < 0, x i+ i > 0)(20) 

|^z+l| (xi — 0, 0) 

0 (x i+ i = 0) 


where 


a = 


K(0+) 


A(0+) + A(0“) ' 

Then, SDEs dTJ) and © can be rewritten as 

dK, 


(13) 


Here we use the fact 0 < a < 1 and equations (HI and 
(usd. Now the meaning of A L becomes clear. Since AH < 
0 in our case (see section E©), one can find A L > 0 when 
XiXi+i < 0, and A L = 0 when XiXi+i > 0. Therefore, 
it is confirmed that part icles g ain ene rgy when they pass 
through the shock front (Zhang, 2000j l. 


dy = s(x) 


>(x) 


dx 


dt + V2KdW 


(14) 


du = - ( + Psynl\ dt ~^^\- dX ~ S Hy) d V\ '( 15 ) 


Derivation of equa tions (THl) and (fj~5l) are the same way 
as of Zhana ( 2000 1. These equations do not contain 8- 
functions and can be integrated directly. Once y(t) is ob¬ 
tained, the position of electrons x(t) can be obtained by 


1/a (y < 0) 

x = ys~\y)=yx{ 2 (y = 0) 

1/(1-a) (y> 0) 


(16) 


In order to see the effect of diffusion, the spatial step size 
of diffusion in one time step At must be larger than that 
of convection, that is vAt < y/2KAt. Hence we derive the 
requirement of the time step as 


At < 


2 A 


(17) 


Functions x(t) and u(t) are numerically integrated as 
follows. We define = X(ti) and X i+ \ = X(U + At), 
where X = x, y and u. First, we discretize equation m 
as 


2/i+l = Vi(Xi) + s(Xi) 


v ( x i ) + ) At 


dx 


+ ^/2K(x hUl )AW 


(18) 


where AW is independent and identically distributed nor¬ 
mal random variable with expected value zero and variance 
At. Then, yi+\ is obtained for given Xi and u,;. Second, 
we get Xi+i from equation (USD- Finally, Ui +1 is calculated 
from equation (USD, which is discretized as 


Ui-\-\ — Ui - 


1 dv c 


3 dx 


(Xi) + p s y n "i{ui) ) At + AT ,(19) 


where 


A L = - 


AH 


\{x i+ i - x^ - s i {yi)(y i +i - yi)\ 


3A K 
AH 


3A Ii 


[xi+i — s 1 {y i )y i+ 1 ] 


-AH 


2.3. Particle Splitting 

A particle splitting is necessary to achieve a wi de mo- _ 

mentu m range of accelerated electrons. Followine lYoshida fc Yanaa 
(il994h . we set splitting surfaces u n in momentum space 
(u s o < u < u s i) with an equal spacing in logarithmic scale: 


u n = u s0 + nAu 


( 21 ) 


where n = 0,1, 2, • • •, n max and A u = (u sl - 'u s0 )/n max . 
Each time an accelerated particle hits the surface u n . the 
particle is split into w particles with the same energy and 
spatial position which particles have attained. The statis¬ 
tical weight which is needed to calculate the final spectrum 
of the particles is decreased by a factor of w in each split¬ 
ting. 


3. Simulations for stationary shock cases 


3.1. Simulation setup 

In the following, we consider electron acceleration at 
stationary plane parallel shock, that is, 


, s Vl+V 2 

Vc( x ) = 2 

(22) 

and AH = v 2 — Ui, so that 


v(x) = l Vl {X< ° ] 

\ V 2 (x > 0) ’ 

(23) 

where constants v\ and v 2 are upstream and downstream 
velocities, respectively. Compression ratio r = vi/v 2 is 
fixed to be 4 throughout the paper. We also assume that 
the diffusion coefficient is uniform both in upstream and 
downstream regions, that is, 

Kc(x) = I Sl±^A , 

(24) 

and A A = I \ 2 — K\, so that 


o' o' 

V A 

II 

(25) 


where upstream and downstream coefficients, K\ and K%, 
depend on electron momentum. In this paper, we assume 


A'i(p) = rK 2 (p) 


= 1.6 x 10 19 B~q 


cm 2 s 1 , (26) 


3[A'(0+) + A'(0 - )] 
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where is the magnetic field strength in units of /iG. 

We set a free escape boundary at x = — Xf e b (< 0) in the 
upstream region. Once a particle goes beyond the bound¬ 
ary, it never comes back to the acceleration site. This fact 
becomes significant when the particle’s penetration depth 
in the upstream region, K\[p)/vi, is comparable to Xf e b- 
On the other hand, when Sf e b is sufficiently large (i.e., 
irfeb —> oo), the particle escape does not occur. 

The start and end times of electron acceleration are 


t a cc = tcooi, where f coo i is the synchrotron cooling time 
given by 


tcool — Psyn ( 


p 

m e c 


-l 


(30) 


and m, we derive 


Using equations 

Pm,cool = ( 6.05 x 10 lb B~Qv‘l)l^m e c . ( 31 ) 

It may happen that th e maxi mum energy is limited 


t = 0 and t ase , respectively. During this period, electrons , ' n,, . , nn I 

• • , i r n by the escape process (Ohira et al., 2010). Characteristic 


are injected with a momentum pj n j at the shock front x = 0 
continuously at a constant rate. Taking an ensemble av¬ 
erage over a number of realizations of SDEs, (HI, ED 
and ED< we obtain the momentum spectrum of the whole 
region (including both upstream and downstream regions 
as well as the shock front) at t = t age . In the following, 
we consider the case (3 > 0. Then, K{p) increases with 
p. Hence, if the injection momentum pj n j is taken so as 
to At < 2K(pi n j)/vf, the condition for time step, equa¬ 
tion urn is always satisfied. It can be seen from Table 2 
that pi n j satisfies this requirement. 

We perform simulations for various parameter sets. Adopted 
parameters are summarized in Tables 1 and 2. The cutoff 
shape of electron spectrum depends on how the maximum 
momentum of electrons is determined. In the next sec¬ 
tion, we consider three cases, age-limited, cooling-limited, 
and escape-limited cases, in order to decide the maximum 
attainable electron momentum due to the diffusive shock 
acceleration. 


spatial length of particles penetrating into the upstream 
region is given by Ki[p)/v\. As long as K\{p)/v i <C atfeb, 
the particles are confined without the significant escape 
loss, and they are accelerated to higher energies. On the 
other hand, when their momentum increases up to suf¬ 
ficiently high energies satisfying K\{p)/v\ > Xf e b, their 
acceleration ceases and they escape into the far upstream. 
Therefore, the maximum momentum of accelerated parti¬ 
cles in this scenario is given by the condition 


Ki{p) 


Vi 

This reads 


— *£feb 


Pm,esc = ( 6.25 X 10 3 B^ 0 X 5 X 15 ) ? m e c 


(32) 


(33) 


where X 15 = Xf e b/10 15 cm. 


3.2. Estimate of Maximum electron momentum 

The maximum momentum of accelerated electrons is 
limited by a fini t e sho c k age, their c ooling or escape (e.g., 


When p m , a ge is smallest among p m , age , p m ,cooi andp miesc , 
the acceleration is limited by the age of the shock. In the 
cooling-limited and escape-limited cases, p m , CO oi and p m , eS c 
are smallest, respectively. 


lYamazaki et al. . 2006p Ohira et al. . 201211 . It is obtained 


by comparisons of timescales, which are given as functions 


3.3. Results 

Figures mm and [3] show results of numerical simula¬ 
tion for the age-limited, cooling-limited and escape-limited 
cases, respectively. The spectra in these figures are for all 


of electron momentum and the shock age, t age in the age- particles which are still in the system at t age , that is 
limited and cooling-limited cases. The acceleration time of 
the diffusive shock acceleration is represented by ([Drury, 

19831 ) 


F{p) oc p 3 /(p) oc 


F(X, It, t ag e) dx 


(34) 


t 


acc 


3 (Ki + K 2\ 
Vl — V2 V Vl V2 J 


(27) 


Using V 2 = V\/r and equation ED with r = 4, we obtain 


= 1.28 X io 4 h ;^ 8 - 2 



( 28 ) 


where vg is the shock velocity, v\, in units of 10 8 cm s _1 . 
From the condition f acc = f age , the age-limited maximum 
momentum, p m , ag e> is derived as 

Pm,age = (2.46 x 10 5 B^cvlt 100 ) Pm e c , (29) 


where fi 00 = t age /100 yr. 

When the magnetic field is strong, the electron accel¬ 
eration is limited by synchrotron cooling. We obtain the 
cooling-limited maximum momentum from the condition 


The value F of the distribution function for given momen¬ 
tum range [it, u + Art] is derived by an ensemble average. 
If injected particles with momentum pi n j have a statisti¬ 
cal weight of unity, then, particles which have experienced 
splitting n times have the statistical weight w~ n . Hence, 
we obtain 

F= Y w ~ n . (35) 

particles 

where summation is taken for all (split) particles which 
have a momentum between u and u + Ait. Error bars 
in spectra of figures 1 -3 are calculated assuming Poisson 
statistics. Taking into account the propagation of errors, 
the statistical error A F for given momentum range [it, u + 
Ait] is calculated as 

(A Ff = Y w ~ 2n ■ ( 36 ) 

particles 
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Table 1: Adopted parameters in the present study. 


Run^ 

P 

B 

[a*g] 

Vl 

[10 8 cm s -1 ] 

^age 

[yr] 

*£feb 

[10 15 cm] 

Pinj 

[m e c] 

Pm,ag^ 

[m e c] 

Pm,coo^ 

[: m e c] 

Pm,es<^ 

[m e c] 

A07-1 

0.7 

1 

1 

3 

oo 

10 3 

3.4 x 10 5 

7.4 x 10 9 

OO 

A07-2 

0.7 

1 

1 

10 

oo 

10 3 

1.9 x 10 6 

7.4 x 10 9 

oo 

A07-3 

0.7 

1 

1 

30 

oo 

10 3 

9.0 x 10 6 

7.4 x 10 9 

oo 

A07-4 

0.7 

1 

1 

100 

oo 

10 3 

5.0 x 10 7 

7.4 x 10 9 

oo 

A07-5 

0.7 

1 

1 

300 

oo 

10 3 

2.4 x 10 8 

7.4 x 10 9 

oo 

A10-1 

1.0 

1 

6 

3 

oo 

10 3 

2.7 x 10 5 

1.5 x 10 9 

oo 

A10-2 

1.0 

1 

6 

10 

oo 

10 3 

8.9 x 10 5 

1.5 x 10 9 

oo 

A10-3 

1.0 

1 

6 

30 

oo 

10 3 

2.7 x 10 6 

1.5 x 10 9 

oo 

A10-4 

1.0 

1 

6 

100 

oo 

10 3 

8.9 x 10 6 

1.5 x 10 9 

oo 

A10-5 

1.0 

1 

6 

300 

oo 

10 3 

2.7 x 10 7 

1.5 x 10 9 

oo 

A15-1 

1.5 

5 

8 

3 

oo 

10 3 

1.8 x 10 4 

1.4 x 10 7 

oo 

A15-2 

1.5 

5 

8 

10 

oo 

10 3 

4.0 x 10 4 

1.4 x 10 7 

oo 

A15-3 

1.5 

5 

8 

30 

oo 

10 3 

8.2 x 10 4 

1.4 x 10 7 

oo 

A15-4 

1.5 

5 

8 

100 

oo 

10 3 

1.8 x 10 5 

1.4 x 10 7 

oo 

A15-5 

1.5 

5 

8 

300 

oo 

10 3 

3.8 x 10 5 

1.4 x 10 7 

oo 

C07-1 

0.7 

2000 

0.1 

10 

oo 

10 4 

1.4 x 10 8 

5.7 x 10 B 

oo 

C07-2 

0.7 

500 

0.1 

100 

oo 

10 4 

5.0 x 10 8 

1.3 x 10 7 

oo 

C10-1 

1.0 

2000 

1 

10 

oo 

10 4 

4.9 x 10 7 

5.5 x 10 6 

oo 

C10-2 

1.0 

500 

1 

100 

oo 

10 4 

1.2 x 10 s 

1.1 x 10 7 

oo 

C15-1 

1.5 

2000 

10 

100 

oo 

10 4 

1.3 x 10 7 

1.6 x 10 B 

oo 

C15-2 

1.5 

500 

10 

800 

oo 

10 4 

2.1 x 10 7 

2.7 x 10 6 

oo 

E07-1 

0.7 

1 

6 

95 

0.1 

10 2 

7.8 x 10 y 

6.1 x 10 1U 

1.3 x 10 5 

E07-2 

0.7 

1 

6 

95 

1 

10 2 

7.8 x 10 9 

6.1 x 10 10 

3.4 x 10 6 

E07-3 

0.7 

1 

6 

95 

10 

10 2 

7.8 x 10 9 

6.1 x 10 10 

9.2 x 10 7 

E07-4 

0.7 

1 

6 

95 

100 

10 2 

7.8 x 10 9 

6.1 x 10 10 

2.5 x 10 9 

E10-1 

1.0 

1 

6 

95 

0.1 

10 2 

8.4 x 10 6 

1.5 x 10 9 

3.8 x 10 3 

E10-2 

1.0 

1 

6 

95 

1 

10 2 

8.4 x 10 6 

1.5 x 10 9 

3.8 x 10 4 

E10-3 

1.0 

1 

6 

95 

10 

10 2 

8.4 x 10 6 

1.5 x 10 9 

3.8 x 10 5 

El 0-4 

1.0 

1 

6 

95 

100 

10 2 

8.4 x 10 6 

1.5 x 10 9 

3.8 x 10 6 

E15-1 

1.5 

1 

6 

95 

0.1 

10 2 

4.1 x 10 4 

2.2 x 10 7 

2.4 x 10 2 

E15-2 

1.5 

1 

6 

95 

l 

10 2 

4.1 x 10 4 

2.2 x 10 7 

1.1 x 10 3 

E15-3 

1.5 

1 

6 

95 

10 

10 2 

4.1 x 10 4 

2.2 x 10 7 

5.2 x 10 3 

E15-4 

1.5 

1 

6 

95 

100 

10 2 

4.1 x 10 4 

2.2 x 10 7 

2.4 x 10 4 


a A, C and E stand for age, cooling and escape, respectively, 
^calculated according to equation ( 1291 ) . 

Calculated according to equation (131I I. 

^calculated according to equation 03j. 
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We adopt different values of u s i, n max and w for different 
runs so as to obtain good statistics near the maximum mo¬ 
mentum (see Appendix). Therefore, for different runs, the 
length of error bars is different from each other. For each 
case, the dependence of the cutoff shape of the electron 
spectrum on j3 is discussed below. 


3.3.1. Age-limited case: p m , age < min{p miCOol ,p m)esc } 

For each value of /? (0.7, 1.0 and 1.5), we had five 
runs with different t age . One can see from figure [1] that at 
lower electron momentum where p -C p m ,age, the spectrum 
is well described by the analytical solution in the steady 
state, F(p) = p 3 /(p) oc p _1 . Furthermore, estimated val¬ 
ues of Pm,age (see Table 1) using equation (1291) agree with 
the simulation result. Hence our present numerical scheme 
works well. 

For all runs of age-limited acceleration, we have fitted 
the spectrum by the following function 


F(p) oc p 1 exp 



(37) 


and obtain the value of cutoff shape parameter a. The 
left panel of figure 0] shows the result. One can roughly 
confirm earlier result of numerical s i mulat ions, a ~ 2/3 
( Kato fe Takaharal l20Q.lt Kang et al. . 2009h . 


3.3.2. Cooling-limited case: p m ,cooi < min{p mag e,p mi esc} 

For each value of (3 (0.7, 1.0 and 1.5), we performed two 
simulations with the magnetic field strength of B = 2 mG 
and 0.5 mG. Since the magnetic field strength B is much 
larger than any other cases, the requirement for time step, 
equation (ED, is the most severe, so that we set larger pi n j 
of 10 4 ?n e c in order to save the computation time. 

In all runs, one can identify the cooling break, at which 
the spectral slope changes from p ~ 1 to p~ 2 . The break en¬ 
ergy p b is determined by the condition t coo i ss i age (jLonsair , 
20111 ) . and we derive 


Pb ~ 2.45 x 10 11 H /i Qt lo 1 o TO e c . 


(38) 


For example, we obtain from this equation pb = 6.1 x 
10 5 m e c for run C07-1, which is consistent with the simu¬ 
lation result. In some runs such as C15-1 and C15-2, we 
can see pile-ups at the high-energy end (see the right panel 
of figured- Based on these facts, we fit spectra which were 
derived from simulations by the following function 


F(p) oc p 1 C b (p) C p (p) exp 


where 


C'b(p) = 1 

and 


C P (P) = 


1 





(39) 


(40) 


(41) 


describe the cooling break and the pile-up effect, respec¬ 
tively. Middle panel of figure @] shows the fitted a as 
a function of f3. The result is roughly consistent with 
the analytical result, a =3 + 1, w hich is derived on the 
steady state assumpti on (IZirakashvili fc AharonianL 12007: 
Yamazaki et all 1201.31) . In the case of (3 = 1.5, fitted value 
of a deviates from the analytical expectation of 2.5. This 
comes from the appearance of pile-up, which deforms the 
spectrum around the maximum momentum. Hence, it is 
implied that the analytical expectation a = (3 + 1 is not 
always hold when the electron pile-up becomes significant. 

P reviously, based on si milar numerical simulation to 
ours, iMarcowith fe Casse ( 201C)h numerically obtain the 
shock front energy spectra, F 0 (p) = F(x = 0,p), for the 
cases of /? = 0, 1/2 and 1, and confirmed the relation 
a = (3 + 1 (for smaller dynamic range of momentum than 
ours). In the present study, we find that at least f3 < 1, 
the relation a = /3 + 1 roughly holds even for the spectrum 
of the whole region. 


3.3.3. Escape-limited case: p m , es c < nhn{p m!age ,p m)COO i} 
In this case, steady-state spectrum at t he shock front 


x = 0 ) has been analytically derived as (Canrioli et al 


2009lj : iReville et all l2009p lYamazaki et al 


2013 1 


F 0 (p) = F(x = 0,p) 


oc p exp 


4 ry(p) (-[ i 0 g y 


1 - e~ l ! y 


(42) 


where y(jp) = (p/pm)^• For each value of (3 (0.7, 1.0 and 
1.5), we have four runs with different xt e b- Then, we find 
that in all runs, the derived spectra are well fitted with 
models described by equation C71) (see figure 0 ■ Hence, 
equation (1421) well reproduces the spectrum of the whole 
region with good accuracy although it has been derived as 
the shock front spectrum at x = 0 . 

Here we discuss on whether the simulated spectra are 
fitted with phenomenological formula, equation (1371) . In 
the limit p p m , we can approximate 1 — e~ 1//y ~ 1 /p 
in equation (H21) . resulting in F n (p) oc exp[— (p/p m )^ ], so 
that one might expect a = f3 ( Yamazaki et al. . 2013 ’!. If 
we fit numerically derived spectra with equation ED, then 
the fitted values of a are not along with the expectation 
a = (3. Difference between equations © and ED with 
a = f3 is large at p ~ p m . However, they seem to lie ona« 
/3 + 0.5 (see the right panel of figure [ID- This implies that 
F(p) oc p _1 exp[— (p/pm)^] is not a good approximation 
around the maximum momentum p m . 


4. Summary and Discussion 

We have developed a numerical method of SDE to sim¬ 
ulate electron acceleration at astrophysical shocks. Our 
code involves Zhang’s method of skew Brownian motion 
and particle splitting. Using this code, we have performed 
simulations of electron acceleration at stationary plane 
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Figure 1: Electron spectra in the age-limited cases with (a) /3 = 0.7 (Runs A07-1, A07-2, A07-3, A07-4, A07-5 from left to right), (b) (3 = 1.0 
(Runs A10-1, A10-2, A10-3, A10-4, A10-5 from left to right) and (c) (3 = 1.5 (Runs A15-1, A15-2, A15-3, A15-4, A15-5 from left to right). 
Lines indicate the best fitted models described by equation 03 - 





Figure 2: Electron spectra in the cooling-limited cases with (a) fi = 0.7 (Runs C07-1, C07-2 from left to right), (b) /3 = 1.0 (Runs C10-1, C10- 
2 from left to right) and (c) /3 = 1.5 (Runs C15-1, C15-2 from left to right). Lines indicate the best fitted models described by equations (1391) . 
(1101) and I-111) . 
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Figure 3: Electron spectra in the escape-limited cases with (a) /3 = 0.7 (Runs E07-1, E07-2, E07-3, E07-4 from left to right), (b) /3 = 1.0 
(Runs E10-1, E10-2, E10-3, E10-4 from left to right) and (c) f3 = 1.5 (Runs E15-1, E15-2, E15-3, E15-4 from left to right). Lines indicate the 
best fitted models described by equation GIJ. 





beta beta beta 


Figure 4: Cutoff shape parameter a as a function of f3 in the age-limited (left), cooling-limited (center) and escape-limited (right) cases. The 
values of a are derived by fitting the simulated spectra with the model described by equations no and for age-limited and cooling-limited 
cases, respectively. Note that for escape-limited cases, we use phenomenological formula, equation J37l) . to derive values of a, while we used 
in Figure 3 the analytical stationary solution, equation to fit the simulated spectra. In the left panel, data points are artificially shifted 
a little in the horizontal direction in order to be separated with each other and to be seen clearly. Solid lines represent a = 2/5, a = (3 +1 and 
a = (3 for the age-limited, cooling-limited and escape-limited cases, respectively. The dashed line in the right panel shows a = [3 + 0.5. 














parallel shock, and we reproduced the analytical result in 
the momentum range much smaller than the maximum 
momentum — f(p ) oc p~ 4 in the age-limited and escape- 
limited cases, and the broken power-law which changes 
from f{p) oc p~ 4 to p~ 5 at the cooling break in the cooling- 
limited cases. These results can be achieved due to incor¬ 
poration of Zhang’s method. Furthermore, the maximum 
electron momentum in the simulated spectra can be well 
explained by simple analytical argument, which is the out¬ 
come of the particle splitting method. Therefore, we be¬ 
lieve our numerical code works well, and it enables us to 
study the cutoff shape of the electron spectrum. 

We have performed simulations for various parameter 
sets, and studied how the cutoff shape, which is charac¬ 
terized by cutoff shape parameter a, changes with the mo¬ 
mentum dependence of the diffusion coefficient /3. In the 
age-limited cases, we have reproduced previous results of 
other authors, a ss 2/3. In the cooling-limited cases, the 
analytical expectation a « f3 + 1 is roughly reproduced al¬ 
though we recognize deviations to some extent (runs Cl5-1 
and C15-2) when the pile-up effect is significant. How¬ 
ever, we have found in the Bohm type diffusion, K oc p 
and (0=1, the cutoff shape parameter a is consistent 
with the analytical prediction a = 2.0 both in the age- 
limited and cooling-limited cases. Hence, if the effect of 
escape can be neglected, a = 2 should be canonical value. 
Note that in the present study, we have assumed plane 
shock geometry and constant electron injection. In real¬ 
ity, the SNR s hock is near ly spherical although it has fluc¬ 
tuation (e.g., Inoue et ah, |2012). In the spherical shock 
case, accelerated particles downstream of the shock expe¬ 
rience adiabatic losses. This reduces the mean energy gain 
they experi ence at the shock, which steepens the spe ctral 
slope (e.g., Yamazaki et, al. . 200(J : Schure et all 201fllh In 
addition, the exact spectral slope in a time-dependent cal¬ 
culation depends on the injection history, which may be 
complicated depending on the shock velocity, ambient den¬ 
sity and so on. These effects may influence the electron 
spectrum. However, we are currently interested in the en¬ 
ergy region near the upper end of the spectrum. At the 
given epoch, the spectrum around the maximum energy 
is dominated by those which are being accelerated at that 
time because previously accelerated particles have suffered 
the adiabatic losses during transported downstream of the 
shock. Hence, one can expect that the cutoff shape of the 
spectrum do es not so much depend on th e past accelera¬ 
tion history ( Yamazaki et al. . 20061 2013h . This issue has 
not yet been studied in detail except for a few works in 
which the spherical and planar sh ock cases were compared 
( Schure et al. . 201(1 Kane . 2015), and should be investi¬ 
gated more in future works. 

The maximum momentum is sometimes determined by 
the escape of accelerated particles upstream. In this case, 
we should use the functional form given by equation flUD , 
otherwise we should use © with a = (3 + 0.5. Electron 
acceleration at SNRs is someti mes limited by the escape 
(see figure 1 and 2 of Ohira et al. . 2012h . as well as proton 


acceleration, which might be i nferred by recent gamma-ray 
observation (e.g., Ohira et al. . 2010l) . In the present study 


we adopt weak magnetic field in all runs of the escape- 
limited cases, so that synchrotron cooling effect can be 
neglected. Hence, our result for the escape-limited cases 
is applicable to the proton acceleration. The cutoff shape 
around the maximum proton momentum may be studied 
by the precise gamma-ray spectrum which will be taken in 
the near future. 

In the present study, we have used the test particle ap¬ 
proximation, neglecting feedback of the accelerated par¬ 
ticles on to the plasma forming background shock struc¬ 
ture. Not electrons but protons deform the background 
plasma, because they are coupled with each other through 
the waves excited by accelerated protons themselves. Var¬ 
ious authors focus on the feedback processes of acceler- 


ties around the shock, (e.g., Berezhko & Ellison 

1999; 

Maikov & Drurvl. 2001: Kang & Jones, 2005t Vladimirov et all 

2006 

Terasawa et all 2007: Caorioli et all 2009a 

: Yamazaki et al. 

2009 

Zirakashvili & Ptuskin, 2012| Bykov et all 

2014). In 


such cosmic-ray modified shocks, electron acceleration is 
also affected by the shock deformation, and the results may 
be different from those in the test particle limit. These 
studies are remained as a future work. 
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Appendix: Parameters for numerical simulation 

Parameters for numerical simulation are summarized in 
Table 2. Time step, At, must satisfy the condition (fT71) . In 
the present study, it is taken to be smaller than 2 K\ (pi n j) /v\. 

We need four parameters, u s o, u s i, n max and w , to 
carry out particle splitting. In the present study, we set 
it s o = ln(pi„j /m e c). The other parameters are taken dif¬ 
ferently for each run. 
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